﻿Imports System
Imports System.IO
Imports System.IO.Ports
Imports System.Linq
Imports System.Collections.Generic
Imports System.ComponentModel
Imports System.Data
Imports System.Drawing
Imports System.Data.SqlServerCe
Imports System.Text

Namespace EfiBase

    Public Class UTM

        Private Const k0 As Double = 0.9996
        Private cooLat As New GpsCoordenadas.Coordinate(0, GpsCoordenadas.CoordinateType.Latitude)
        Private cooLong As New GpsCoordenadas.Coordinate(0, GpsCoordenadas.CoordinateType.Longitude)
        'Private utmDat As GpsCoordenadas.UTMDatum = GpsCoordenadas.UTMDatum.SIRGAS2000
        'Private utmDat As GpsCoordenadas.UTMDatum = GpsCoordenadas.UTMDatum.SAD69
        'Private utmDat As GPS.UTMDatum = GPS.UTMDatum.WGS84
        Private utmDat As GpsCoordenadas.UTMDatum = GpsCoordenadas.UTMDatum.International1924
        Private dblNorthing As Double = 0
        Private dblEasting As Double = 0
        Private utmL As GpsCoordenadas.UTMLat = GpsCoordenadas.UTMLat.South
        Private intZone As Int32

        Public Property UTMLatitudeHemisphere() As GpsCoordenadas.UTMLat
            Get
                Return utmL
            End Get
            Set(ByVal value As GpsCoordenadas.UTMLat)
                utmL = value
                GetLatLong()
            End Set
        End Property

        Public Property Datum() As GpsCoordenadas.UTMDatum
            Get
                Return utmDat
            End Get
            Set(ByVal value As GpsCoordenadas.UTMDatum)
                utmDat = value
            End Set
        End Property

        Public Property Latitude() As GpsCoordenadas.Coordinate
            Get
                Return cooLat
            End Get
            Set(ByVal value As GpsCoordenadas.Coordinate)
                If value.Direction = GpsCoordenadas.CardinalDirection.North OrElse value.Direction = GpsCoordenadas.CardinalDirection.South Then
                    cooLat = value
                    utmL = DirectCast(cooLat.Direction, GpsCoordenadas.UTMLat)
                    GetUTM()
                End If
            End Set
        End Property

        Public Property Longitude() As GpsCoordenadas.Coordinate
            Get
                Return cooLong
            End Get
            Set(ByVal value As GpsCoordenadas.Coordinate)
                If value.Direction = GpsCoordenadas.CardinalDirection.East OrElse value.Direction = GpsCoordenadas.CardinalDirection.West Then
                    cooLong = value
                    GetUTM()
                End If
            End Set
        End Property

        Public Property Zone() As Int32
            Get
                Return intZone
            End Get
            Set(ByVal value As Int32)
                intZone = value
                GetLatLong()
            End Set
        End Property

        Public Property Easting() As Double
            Get
                Return dblEasting
            End Get
            Set(ByVal value As Double)
                dblEasting = value
                GetLatLong()
            End Set
        End Property

        Public Property Northing() As Double
            Get
                Return dblNorthing
            End Get
            Set(ByVal value As Double)
                dblNorthing = value
                GetLatLong()
            End Set
        End Property

        Public Sub New(ByVal iLatitude As GpsCoordenadas.Coordinate, ByVal iLongitude As GpsCoordenadas.Coordinate)
            cooLat = iLatitude
            cooLong = iLongitude
            intZone = GetZone()
            utmL = DirectCast(cooLat.Direction, GpsCoordenadas.UTMLat)
            GetUTM()
        End Sub

        Public Sub New(ByVal Northing As Double, ByVal Easting As Double, ByVal Zone As Int32)
            dblNorthing = Northing
            dblEasting = Easting
            intZone = Zone
            GetLatLong()
        End Sub

        Public Function GetZone() As Int32
            Dim decLongAbs As Decimal = cooLong.GetDecimalCoordinate
            If cooLong.Direction = GpsCoordenadas.CardinalDirection.West Then
                Return Convert.ToInt32(Math.Ceiling((180 + decLongAbs) / 6))
            Else
                Return Convert.ToInt32(Math.Ceiling(decLongAbs / 6) + 30)
            End If

        End Function

        Public Function GetZoneCM() As Int32
            Return 6 * intZone - 183
        End Function

        Public Sub GetUTM()
            Dim a As Double = LookupA()
            Dim b As Double = LookupB()
            Dim f As Double = (a - b) / a
            Dim invf As Double = 1 / f
            Dim rm As Double = (a * b) ^ 0.5
            Dim e As Double = Math.Sqrt(1 - (b / a) ^ 2)
            Dim e1sq As Double = e ^ 2 / (1 - e ^ 2)
            Dim n As Double = (a - b) / (a + b)
            Dim latRad As Double = cooLat.GetDecimalCoordinate * Math.PI / 180
            Dim rho As Double = a * (1 - e ^ 2) / ((1 - (e * Math.Sin(latRad)) ^ 2) ^ (3 / 2))
            Dim nu As Double = a / ((1 - (e * Math.Sin(latRad)) ^ 2) ^ (1 / 2))

            Dim a0 As Double = a * (1 - n + (5 * n * n / 4) * (1 - n) + (81 * n ^ 4 / 64) * (1 - n))
            Dim b0 As Double = (3 * a * n / 2) * (1 - n - (7 * n * n / 8) * (1 - n) + 55 * n ^ 4 / 64)
            Dim c0 As Double = (15 * a * n * n / 16) * (1 - n + (3 * n * n / 4) * (1 - n))
            Dim d0 As Double = (35 * a * n ^ 3 / 48) * (1 - n + 11 * n * n / 16)
            Dim e0 As Double = (315 * a * n ^ 4 / 51) * (1 - n)
            Dim s As Double = a0 * latRad - b0 * Math.Sin(2 * latRad) + c0 * Math.Sin(4 * latRad) - d0 * Math.Sin(6 * latRad) + e0 * Math.Sin(8 * latRad)

            Dim p As Double = (cooLong.GetDecimalCoordinate - GetZoneCM()) * 3600 / 10000
            Dim sin1 As Double = Math.PI / (180 * 3600)

            Dim ki As Double = s * k0
            Dim kii As Double = nu * Math.Sin(latRad) * Math.Cos(latRad) * sin1 ^ 2 * k0 * (100000000) / 2
            Dim kiii As Double = ((sin1 ^ 4 * nu * Math.Sin(latRad) * Math.Cos(latRad) ^ 3) / 24) * (5 - Math.Tan(latRad) ^ 2 + 9 * e1sq * Math.Cos(latRad) ^ 2 + 4 * e1sq ^ 2 * Math.Cos(latRad) ^ 4) * k0 * (10000000000000000)
            Dim kiv As Double = nu * Math.Cos(latRad) * sin1 * k0 * 10000
            Dim kv As Double = (sin1 * Math.Cos(latRad)) ^ 3 * (nu / 6) * (1 - Math.Tan(latRad) ^ 2 + e1sq * Math.Cos(latRad) ^ 2) * k0 * (1000000000000)
            Dim a6 As Double = ((p * sin1) ^ 6 * nu * Math.Sin(latRad) * Math.Cos(latRad) ^ 5 / 720) * (61 - 58 * Math.Tan(latRad) ^ 2 + Math.Tan(latRad) ^ 4 + 270 * e1sq * Math.Cos(latRad) ^ 2 - 330 * e1sq * Math.Sin(latRad) ^ 2) * k0 * (1.0E+24)

            dblEasting = 500000 + (kiv * p + kv * p ^ 3)
            dblNorthing = (ki + kii * p * p + kiii * p ^ 4)
            If cooLat.Direction = GpsCoordenadas.CardinalDirection.South Then dblNorthing = 10000000 + dblNorthing
        End Sub

        Public Sub GetLatLong()
            Dim a As Double = LookupA()
            Dim b As Double = LookupB()
            Dim f As Double = (a - b) / a
            Dim invf As Double = 1 / f
            Dim rm As Double = (a * b) ^ 0.5
            Dim ec As Double = Math.Sqrt(1 - (b / a) ^ 2)
            Dim eisq As Double = ec ^ 2 / (1 - ec ^ 2)

            dblNorthing = If(utmL = GpsCoordenadas.UTMLat.North, dblNorthing, 10000000 - dblNorthing)

            Dim arc As Double = dblNorthing / k0
            Dim mu As Double = arc / (a * (1 - ec ^ 2 / 4 - 3 * ec ^ 4 / 64 - 5 * ec ^ 6 / 256))
            Dim ei As Double = (1 - (1 - ec * ec) ^ (1 / 2)) / (1 + (1 - ec * ec) ^ (1 / 2))
            Dim ca As Double = 3 * ei / 2 - 27 * ei ^ 3 / 32
            Dim cb As Double = 21 * ei ^ 2 / 16 - 55 * ei ^ 4 / 32
            Dim ccc As Double = 151 * ei ^ 3 / 96
            Dim cd As Double = 1097 * ei ^ 4 / 512
            Dim phi1 As Double = mu + ca * Math.Sin(2 * mu) + cb * Math.Sin(4 * mu) + ccc * Math.Sin(6 * mu) + cd * Math.Sin(8 * mu)

            Dim sin1 As Double = Math.PI / (180 * 3600)
            Dim q0 As Double = eisq * Math.Cos(phi1) ^ 2
            Dim t0 As Double = Math.Tan(phi1) ^ 2
            Dim n0 As Double = a / (1 - (ec * Math.Sin(phi1)) ^ 2) ^ (1 / 2)
            Dim r0 As Double = a * (1 - ec * ec) / (1 - (ec * Math.Sin(phi1)) ^ 2) ^ (3 / 2)
            Dim dd0 As Double = (500000 - dblEasting) / (n0 * k0)

            Dim fact1 As Double = n0 * Math.Tan(phi1) / r0
            Dim fact2 As Double = dd0 * dd0 / 2
            Dim fact3 As Double = (5 + 3 * t0 + 10 * q0 - 4 * q0 * q0 - 9 * eisq) * dd0 ^ 4 / 24
            Dim fact4 As Double = (61 + 90 * t0 + 298 * q0 + 45 * t0 * t0 - 252 * eisq - 3 * q0 * q0) * dd0 ^ 6 / 720

            Dim lof1 As Double = dd0
            Dim lof2 As Double = (1 + 2 * t0 + q0) * dd0 ^ 3 / 6
            Dim lof3 As Double = (5 - 2 * q0 + 28 * t0 - 3 * q0 ^ 2 + 8 * eisq + 24 * t0 ^ 2) * dd0 ^ 5 / 120

            Dim lat As Double = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / Math.PI
            If utmL = GpsCoordenadas.UTMLat.South Then lat = -lat
            cooLat = New GpsCoordenadas.Coordinate(Convert.ToDecimal(lat), GpsCoordenadas.CoordinateType.Latitude)

            Dim lon As Double = GetZoneCM() - ((lof1 - lof2 + lof3) / Math.Cos(phi1)) * 180 / Math.PI

            cooLong = New GpsCoordenadas.Coordinate(Convert.ToDecimal(lon), GpsCoordenadas.CoordinateType.Longitude)
        End Sub

        Private Function LookupA() As Double
            Select Case utmDat
                Case GpsCoordenadas.UTMDatum.SIRGAS2000
                    Return 6378137
                Case GpsCoordenadas.UTMDatum.SAD69
                    Return 6378160
                Case GpsCoordenadas.UTMDatum.NAD83, GpsCoordenadas.UTMDatum.WGS84, GpsCoordenadas.UTMDatum.GRS80
                    Return 6378137
                Case GpsCoordenadas.UTMDatum.WGS72
                    Return 6378135
                Case GpsCoordenadas.UTMDatum.Australian1965
                    Return 6378160
                Case GpsCoordenadas.UTMDatum.Krasovsky1940
                    Return 6378245
                Case GpsCoordenadas.UTMDatum.NorthAmerican1927
                    Return 6378206.4
                Case GpsCoordenadas.UTMDatum.International1924, GpsCoordenadas.UTMDatum.Hayford1909
                    Return 6378388
                Case GpsCoordenadas.UTMDatum.Clarke1880
                    Return 6378249.1
                Case GpsCoordenadas.UTMDatum.Clarke1866
                    Return 6378206.4
                Case GpsCoordenadas.UTMDatum.Bessel1841
                    Return 6377397.2
                Case GpsCoordenadas.UTMDatum.Airy1830
                    Return 6377563.4
                Case GpsCoordenadas.UTMDatum.Everest1830
                    Return 6377276.3
            End Select
        End Function

        Private Function LookupB() As Double
            Select Case utmDat
                Case GpsCoordenadas.UTMDatum.SIRGAS2000
                    Return 6356752.31414
                Case GpsCoordenadas.UTMDatum.SAD69
                    Return 6356774.7192
                Case GpsCoordenadas.UTMDatum.NAD83, GpsCoordenadas.UTMDatum.WGS84
                    Return 6356752.3142
                Case GpsCoordenadas.UTMDatum.GRS80
                    Return 6356752.31414
                Case GpsCoordenadas.UTMDatum.WGS72
                    Return 6356750.5
                Case GpsCoordenadas.UTMDatum.Australian1965
                    Return 6356774.7
                Case GpsCoordenadas.UTMDatum.Krasovsky1940
                    Return 6356863
                Case GpsCoordenadas.UTMDatum.NorthAmerican1927
                    Return 6356583.8
                Case GpsCoordenadas.UTMDatum.International1924, GpsCoordenadas.UTMDatum.Hayford1909
                    Return 6356911.9
                Case GpsCoordenadas.UTMDatum.Clarke1880
                    Return 6356514.9
                Case GpsCoordenadas.UTMDatum.Clarke1866
                    Return 6356583.8
                Case GpsCoordenadas.UTMDatum.Bessel1841
                    Return 6356079
                Case GpsCoordenadas.UTMDatum.Airy1830
                    Return 6356256.9
                Case GpsCoordenadas.UTMDatum.Everest1830
                    Return 6356075.4
            End Select
        End Function

        Public Function GetZoneChar() As Char
            Select Case cooLat.GetDecimalCoordinate
                Case -90 To -84
                    Return "A"c
                Case -84 To -72
                    Return "C"c
                Case -72 To -64
                    Return "D"c
                Case -64 To -56
                    Return "E"c
                Case -56 To -48
                    Return "F"c
                Case -48 To -40
                    Return "G"c
                Case -40 To -32
                    Return "H"c
                Case -32 To -24
                    Return "J"c
                Case -24 To -16
                    Return "K"c
                Case -16 To -8
                    Return "L"c
                Case -8 To 0
                    Return "M"c
                Case 0 To 8
                    Return "N"c
                Case 8 To 16
                    Return "P"c
                Case 16 To 24
                    Return "Q"c
                Case 24 To 32
                    Return "R"c
                Case 32 To 40
                    Return "S"c
                Case 40 To 48
                    Return "T"c
                Case 48 To 56
                    Return "U"c
                Case 56 To 64
                    Return "V"c
                Case 64 To 72
                    Return "W"c
                Case 72 To 84
                    Return "X"c
                Case Else
                    Return "Z"c
            End Select
        End Function

    End Class

End Namespace